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Abstract. Graphene antidot lattices constitute a novel class of nano-engineered 
graphene devices with controllable electronic and optical properties. An antidot lattice 
consists of a periodic array of holes which causes a band gap to open up around the 
Fermi level, turning graphene from a semimetal into a semiconductor. We calculate 
the electronic band structure of graphene antidot lattices using three numerical 
approaches with different levels of computational complexity, efficiency, and accuracy. 
Fast finite-element solutions of the Dirac equation capture qualitative features of the 
band structure, while full tight-binding calculations and density functional theory are 
necessary for more reliable predictions of the band structure. We compare the three 
computational approaches and investigate the role of hydrogen passivation within our 
density functional theory scheme. 
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1. Introduction 

Since its discovery in 2004 [U [2] , graphene has become a research field of tremendous 
interest within the sohd state physics community [3]. The interest stems from the 
particular electronic properties of graphene as well as the promising perspectives for 
future technological applications [Ij. The electronic excitations around the Fermi level 
of graphene resemble those of massless, relativistic Dirac fermions, allowing predictions 
from quantum electrodynamics to be tested in a solid state system [5]. From a 
technological point of view, several future applications have already been envisioned. 
These include the use of graphene for single molecule gas detection |6], graphene- 
based field-effect transistors [Ij, and quantum information processing in nano-engineered 
graphene sheets [7j. Additionally, graphene is the strongest material ever tested, 
suggesting the use of carbon-fiber reinforcements in novel material composites [8]. 

Metamaterials constitute another popular field of research in contemporary science. 
Contrary to conventional, naturally occurring materials, metamaterials derive their 
properties from their artificial, man-made, periodic small-scale structure rather than 
their chemical or atomic composition [9]. When properly designed and fabricated, 
metamaterials offer optimized and unusual, sometimes even counter-intuitive, responses 
to specific excitations [10]. Examples include metamaterials with negative permittivity 
and permeability [H], superlenses [12], [13], and cloaking devices [H]. Photonic [T5| [T6] 
and phononic [17] crystals are closely related to metamaterials, although they are 
typically designed to alter the response to electromagnetic and acoustic excitations, 
respectively, at wavelengths similar to the dimensions of the small-scale structure. 
The realization of artificial band structures in two-dimensional electron gasses may be 
pursued with similar approaches [TBI [I9], allowing the formation of e.g. Dirac cones in 
conventional antidot lattices [201 HI] • 

Based on the above ideas, some of us have recently proposed to alter in a controllable 
manner the electronic and optical properties of graphene by fabricating a periodic 
arrangement of perforations or holes in a graphene sheet [22]. We refer to this kind of 
structure as a graphene antidot lattice owing to its close resemblance with conventional 
antidot lattices defined on top of a two-dimensional electron gas in a semiconductor 
heterostructure [23l IMj . Using tight-binding calculations we have shown that such 
a periodic array of holes in a graphene sheet causes a band gap to open up around 
the Fermi level ^22], changing graphene from a semimetal to a semiconductor with 
corresponding clear signatures in the optical excitation spectrum [25]. Soon after our 
proposal, graphene antidot lattices were realized experimentally by Shen and co-workers 
[26] and Eroms and Weiss [2^ with lattice constants below 100 nm. The rapidly 
improving ability to pattern monolayer films with e-beam lithography suggests that 
graphene antidot lattices with typical dimensions towards the 10 nm scale may be within 
reach [28l [29] . Furthermore, Girit and co-workers recently monitored the dynamics at 
the edges of a growing hole in real time using a transmission electron microscope [30] . 
and Jia and co-workers demonstrated a method for producing graphitic nanoribbon 
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edges in a controlled manner via Joule heating [21]. Very recently, Rodriguez- Manzo 
and Banhart created individual vacancies in carbon nanotubes using a 1 A diameter 
e-beam [32]. These advances suggest that fabrication of nano-scale graphene antidot 
lattices with desired hole geometries may be possible in the near future. 

In the endeavors of modeling these structures one is faced with a compromise 
between computational efficiency and accuracy. Small-scale lattices with perfect 
periodicity and identical few-nm sized holes can be treated accurately with density 
functional theory (DFT), but this is a computationally heavy and time consuming 
approach, which limits the possibilities to perform large, systematic studies. For 
example, in order to model lattice disorder, such as variations in the hole geometry 
and alignment, it may be necessary to form a super cell containing several holes at the 
cost of an increased computational time. In order to circumvent this problem, one can 
make use of the pseudo-relativistic behavior of electrons in bulk graphene close to the 
Fermi level and solve the corresponding Dirac equation using computationally cheaper 
methods, however, possibly at the cost of a decreased computational accuracy. 

The aim of this paper is to study the band structure of graphene antidot lattices 
using three numerical approaches of different computational complexity, efficiency, and 
accuracy. We first develop a computationally cheap scheme based on a finite-element 
solution of the Dirac equation. This method gives reasonable predictions for the size 
of the band gap due to the antidot lattice, but has limited accuracy in predicting the 
full band structure. For better predictions of the band structure, we employ a vr-orbital 
tight-binding scheme, which is still numerically cheap and capable of treating larger 
antidot lattices. The results are compared with computationally demanding, fuU-fiedged 
ab initio calculations, based on density functional theory, which we expect to predict 
the band structure with the highest accuracy. The tight-binding calculations agree well 
with qualitative features of the band structure calculations based on density functional 
theory, although some differences are found on a quantitative level. Finally, we discuss 
hydrogen passivation along the edges of the holes in a graphene antidot lattice and study 
the infiuence on the electronic properties using density functional theory. 

The paper is organized as follows: In Section [2] we introduce graphene antidot 
lattices and give a brief overview of the existing literature on the topic. In Section [3] 
we describe our three computational approaches; finite-element solutions of the Dirac 
equation (DE), a vr-orbital tight-binding scheme (TB), and density functional theory 
calculations (DFT). A comparison and discussion of the results obtained using the three 
methods are given in Section|H Finally, we discuss in Section[5]the infiuence of hydrogen 
passivation on the band structure, before stating our conclusions in Section [6l 

2. Graphene antidot lattices 

A graphene antidot lattice consists of a periodic arrangement of holes in a graphene 
sheet [22]. In the following, we consider a hexagonal lattice of circular holes, but other 
lattice structures, e.g. square lattices, with different holes shapes are expected to exhibit 
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{12,3} {7,3} {10,6.4} 




Figure 1. Unit cells of three hexagonal graphene antidot lattices with different side 
lengths L and hole radii R. The structures are denoted as {L, R} with both lengths 
measured in units of the graphene lattice constant a ~ 2.46 A. Here we have assumed 
that the edges of the holes have been hydrogen passivated (hydrogen shown as white 
atoms). 



similar physics. In particular, we anticipate an opening of a band gap around the Fermi 
level for a large class of antidot lattices pS] . The hexagonal unit cells with different hole 
sizes are shown in Fig. [H The structures are characterized by the side lengths L of the 
hexagonal unit cells and the approximate radii R of the holes, both measured in units of 
the graphene lattice constant a = y/3lc — 2.46 A, where Ic = 1-42 A is the bond length 
between neighboring carbon atoms. In Fig. [1], the holes are assumed to be passivated 
with hydrogen, using the bond length 1.1 A between neighboring carbon and hydrogen 
atoms. Throughout the paper, we denote a given structure as {L,R}, where L is an 
integer, but R not necessarily. We will consider only very small structures with L ~ 10. 
Although it may not be conceivable to fabricate such small structures within the near 
future, the small unit cells allow for systematic comparisons of our three computational 
schemes. In particular, with small unit cells we can perform computationally heavy 
DFT calculations. Importantly, simple scaling relations have been demonstrated for the 
size of the band gap in terms of the total number of atoms and the number of removed 
atoms within a unit cell, making it possible to extrapolate results to larger geometries 
|22j . Such scaling relations may be helpful when modeling on-going experiments on 
graphene antidot lattices [26| [27]. 

In our original proposal for graphene antidot lattices, we focused on the possibility 
of fabricating intentional 'defects' by leaving out one or more holes in the otherwise 
periodic structure [22]. As we showed, such defects lead to the formation of localized 
electronic states at the locations of the defects with energies inside the band gap. Several 
such (possibly coupled) defects would then form a platform for coupled electronic spin 
qubits in a graphene-based quantum computing architecture [22]. Similar ideas based 
on conventional antidot lattices defined on a two-dimensional gas in a semiconductor 
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heterostructure have previously been studied by some of us [ISl [12] • However, as already 
mentioned, the perfectly periodic graphene antidot lattice constitutes an interesting 
structure on its own. In particular, the controllable opening of a band gap may 
potentially pave the way for graphene-based semiconductor devices. In Ref. [25] some 
of us studied the optical properties of graphene antidot lattices, showing that they 
behave as dip ole- allowed direct gap two-dimensional semiconductors with a pronounced 
optical absorption edge. Additional studies of the electronic properties have been 
performed by Vanevic, Stojanovic, and Kindermann [M] as well as by some of us 
[33] . Vanevic and co-workers studied in detail the occurrence of flat bands due to 
sublattice imbalances and irregularities in the hole shapes at the atomic level. In 
our study, we addressed the roles of geometry relaxation and electron spin using DFT 
calculations. Very recently, Rosales and co-workers studied the transport properties of 
antidot lattices along graphene nanoribbons [35]. Turning around the ideas of making 
graphene semiconducting using periodic superlattices, it has recently been shown that 
periodic potential modulations may create graphene-like electronic band structures of 
two-dimensional gases in semiconductor heterostructures [20l [2T] . In that case, the 
possibility to control the slope of the linear bands and thus the velocity of the Dirac 
fermions is of great interest. 

3. Computational methods 

In the following we outline the three computational methods employed in this work. 
As a computationally cheap approach we consider first finite-element solutions of the 
Dirac equation (DE). Within this approach, large unit cells can be treated and the 
computations are fast. The method relies on the linear bands of bulk graphene 
around the Fermi level. As a more refined approach, we consider next vr-orbital tight- 
binding calculations (TB). This method goes beyond the assumption of a linear band 
structure of bulk graphene, and the edges of the antidot holes can be carefully treated, 
including possible effects due to valley mixing. Finally, we consider full-fledged ab initio 
calculations using DFT. While this method is computationally heavy, DFT is a widely 
used standard for doing first principles calculations and we expect it to provide the most 
detailed description of the electronic structure. 

3.1. Dirac equation (DE) 

We first describe our finite-element solutions of the Dirac equation. The method is based 
on the band structure of bulk graphene close to the two Dirac points being linear and 
well described by the Dirac equation |3]. Within this picture, the atomic honeycomb 
lattice structure of graphene is replaced by an effective continuum description. As an 
example, we show in Figs. [2|i and Wp, respectively, a graphene antidot lattice unit cell 
and the corresponding continuum domain on which the Dirac equation is solved. The 
hole in the unit cell is mimicked with a mass term M(r) in the Dirac equation at the 
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Figure 2. Unit cell, continuum domain, and finite-element mesh, a) Hexagonal unit 
cell of the {7,3} graphene antidot lattice, b) Corresponding continuum domain on 
which the Dirac equation is solved. The hole (hatched area) is modeled with a mass 
term M(r) in the Dirac equation. The normal vector to the hole n, forming the angle 
(j) with the horizontal axis, is used to define appropriate boundary conditions along the 
edge of the hole (see text) c) Corresponding finite-element mesh on which we solve the 
Dirac equation. The edge of the hole is shown with red. Periodic Bloch conditions are 
imposed on the outer boundary of the unit cell. 



location of the hole; see explanation following Eq. ([2]). For large masses, the Dirac 
fermions are effectively excluded from the location of the hole and the mass term can 
be replaced by appropriate boundary conditions along the edge of the hole, indicated 
with red in Fig. [2t. In Fig. ^ we also show an example of the finite-element mesh on 
which the Dirac equation is discretized and solved. Periodic Bloch boundary conditions 
are imposed on the outer edges of the unit cell, making the problem equivalent to that 
of an infinitely large graphene antidot lattice. 

Electronic states close to one of the two Dirac points of bulk graphene can be 
expressed in terms of envelope wave functions contained in the two-component spinor 
|\E') with one component corresponding to each of the two sublattices in the honeycomb 
structure of graphene [3]. Spinors corresponding to states close to one of the Dirac 
points satisfy the Dirac equation 

H\^>) = [VF t> a + M(r)<T,] = E\^>), (1) 

where vp — 10^ ms^"*^ is the Fermi velocity [2J, p = [px^Vy] is the momentum, cr = [a^, d'y] 
is the pseudo-spin corresponding to the two sublattices, and M(r) is the mass that 
couples to dz and is non-zero only inside the holes. Spinors associated with the other 
Dirac point satisfy Eq. ([T]) with the replacement a &* = [ax,— d'y]. Within this 
description, states close to different Dirac points are assumed not to couple. The real- 
space representation of the spinor is \l'(r) = (r|\l') = ['?/;i(r), ■?/'2(i")]"^, where ipi and 
ip2 are the envelope functions corresponding to each of the two sublattices. Equation 
([T]) is correspondingly written 

M(r) —ihvpidx — idy) 

—ihypidx + idy) — M(r) 

We now consider the situation, where Dirac fermions are excluded from the holes, by 
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taking the hmit M(r) oo inside the holes. In that limit, we can derive the appropriate 
boundary conditions for the spinor along the edges of a hole and solve the resulting 
problem outside the holes. The boundary conditions are derived by requiring that no 
particle current runs into a hole. The particle current operator is j = Vp-ff = vpa, 
and the local particle current density in the state "^{r) is j(r) = \E'''"(r)j\l'(r). Imposing 
n ■ j(r) = along the edge of a hole with n being the outward-pointing normal vector 
to the hole, one can derive the condition tpi{r) = ze~*'^^/'2(r) along the boundary, where 
the angle is defined in Fig. [2b. This procedure was originally developed by Berry and 
Mondragon in studies of neutrino billiards [36] and more recently employed by Tworzydlo 
and co-workers in the context of graphene [37j. Along the outer boundaries of the unit 
cell we impose periodic Bloch boundary conditions, and we are thus left with a system 
of coupled differential equations on a finite-size domain with well defined boundary 
conditions. Problems of this type are well suited for commercially available finite- 
element solvers, and the numerical implementation is relatively straightforward and 
fast using the standard finite-element package COMSOL Multiphysics [38]. The finite- 
element solver discretizes and solves the problem on an optimized mesh of the finite-size 
domain. The mesh shown in Fig. [2b was generated with COMSOL Multiphysics. 



3.2. Tight-binding (TB) 

We next describe our tight-binding scheme. The Dirac equation approach introduced 
above is a continuum description of the electronic properties, ignoring the detailed 
atomic structure of graphene and the edges of the holes, which may lead to scattering 
between the two Dirac points. It moreover assumes linear bands of bulk graphene. 
To capture effects of the atomic structure, including the influence of edge geometry, 
and in order to incorporate a realistic description of the band structure of bulk 
graphene, we need to go beyond the simple Dirac fermion picture. In our tight-binding 
scheme, the starting point is the Schrodinger equation for a single electron in real-space 
representation 



-V' + Vir) 



^(r) = e^(r), (3) 



2me 

where V is an effective potential and nie is the electron mass. The unknown eigenstate 
\tp) is subsequently expanded in a set of localized "atomic" wave functions |-R, /) as a 
superposition \ip) = J2^ri\^^^) "with expansion coefficients C^^. Here, each atomic 
state is labeled by the orbital symmetry {l=s,Px,Pz---) and the position of the atom R. 
This transforms the Schrodinger equation into a matrix equation reading 

J2{RJ\H^''\R'J')Cn'M = eJ2{RJ\R'J')Cri>,,- (4) 

R',l' R',V 

At this point, several approximations can be adopted in order to simplify the 
calculations. First, the atomic orbitals are usually taken to be orthogonal, i.e., 
{R,l\R',l') = 6j^j^,6i^i'. This means that the matrix problem becomes a simple rather 
than a generalized eigenvalue problem. Second, the matrix elements of H"^^ are regarded 
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as empirical parameters fitted, usually, to experimental data. In the simplest tight- 
binding description of planar carbon structures contained in the (x, ?/)-plane, just a single 
Pz or TT-orbital on each site is considered and only nearest-neighbor matrix elements are 
retained. This "hopping integral" is denoted as —(3, with j3 ~ 3.033 eV [39]. Other 
values of the hopping integral can also be found in the literature. For example, the choice 
/3 ~ 2.7 eV provides low-energy band structures for bulk graphene consistently with 
density functional theory calculations [^. However, the Fermi velocity is determined 
by the relation vp = \/3Pa/2h and by choosing P ^ 3.033 eV, we obtain vp = 9.9 x 10^ 
ms~^ in good agreement with experiments [2]. 

The reason for considering only vr-orbitals is that vr-orbitals with odd 2;-parity 
decouple from the cr-orbitals spanned by s,Px, and Py states that all have even z-parity. 
Moreover, the bands in the vicinity of the band gap are all produced by the loosely 
bound TT-orbitals. Hence, for all structures considered in the present work, we need only 
include vr-orbitals explicitly. Also, even though realistic structures will contain hydrogen 
terminated edges, the hydrogen atoms couple only to the cx-orbitals and are therefore 
irrelevant for vr-states. In a more sophisticated model, bare or hydrogen terminated 
edges lead to a small modification of the vr-electron hopping integrals near an edge due 
to relaxation of the geometry. This modification is ignored as it simply leads to a small 
additional opening of the band gap [22] . 

3.3. Density functional theory (DFT) 

Finally, we discuss our DFT calculations. This method provides the most detailed 
description of graphene antidot lattices, and we expect it to yield the most accurate 
results. The accuracy comes at the cost of the method being numerically demanding and 
the required computational resources exceed those typically available on a standard PC. 
Density functional theory is a widely used standard for electronic structure calculations 
and we shall here only briefly outline the underlying theory [51]. 

The method takes as starting point the full interacting many-body system involving 
all electrons and atom nuclei making up the graphene antidot lattice. Diagonalizing 
the corresponding many-body Hamiltonian is a tremendous task, but the problem can 
be brought to a somewhat simpler form using the Born-Oppenheimer approximation 
in which the positions of the nuclei are fixed. We are then considering a system of 
interacting electrons moving in an external potential created by the nuclei at fixed 
positions. This is still a difficult many-body problem, but further advances can be 
made following Hohenberg and Kohn who showed that the ground state energy is 
uniquely determined by the ground state electron density [12]. Kohn and Sham (KS) 
later realized that this density can be obtained from a single-particle picture of non- 
interacting electrons. The corresponding Hamiltonian for the single-particle KS orbitals 
ilji is expressed by the KS equations as [13] 



(5) 
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where the effective potential 

depends exphcitly on the density p(r) = P with the sum running over occupied 

KS orbitals. Here, K(r, {Rj^}) is the external potential due to the atoms at positions 
Rj^. The so-called exchange-correlation term Exc{r) accounts for all many-body effects 
and is not known exactly, but must be appropriately approximated. Finally, the ground 
state energy of the interacting problem is 



B[p(r)| = r|p(r)| + j <irp(r)V'<.(r, {R,.}) 



+ Ufdrdrp^ + EJp(r)], (7) 



2 „ „ 

where T is the kinetic energy corresponding to the density p(r). 

We are now left with the problem of determining the density p(r). The density is a 
function of only three coordinates, unlike the A^-particle wavefunction of 3A^ coordinates. 
The set of KS equations is solved self-consistently: starting from an initial density, the 
effective potential is computed together with the KS orbitals and the corresponding 
density, and this procedure is repeated until convergence has been reached. The 
band structure can then by calculated corresponding to the chosen coordinates of the 
nuclei Rj^. The total energy of the system can further be minimized with respect 
to the coordinates of the nuclei. This is referred to as geometry relaxation. The 
method can easily be extended to include spin as well as different species of nuclei. 
In this work, we use spin-polarized DFT as implemented in the Siesta code |44j . 
The structures are relaxed using computationally cheaper DFT based tight-binding 
methods [15]. Performing electronic structure calculations using DFT on geometries 
relaxed in this way is known to provide accurate results [33] . As commonly done, the 
core electrons are replaced by pseudo-potentials and the remaining valence-electrons 
are described with localized atomic orbitals. For the exchange-correlation potential 
we employ the widely used Perdew-Burke-Ernzerhof parametrization of the generalized 
gradient approximation [36]. We mainly use a so-called double-^ polarized (DZP) basis 
set size, consisting of 13 functions per carbon atom. Contrary to the DE and TB 
methods, the antidot edges are hydrogen-passivated in the DFT calculations. The effects 
of passivation are discussed in Section [5l Further details of our DFT calculations can 
be found in Ref. [331. 



4. Band structures 



We now present and compare our results for the electronic band structure obtained using 
the three methods described in the previous section. This provides valuable insight into 
the physics dominating the electronic properties of graphene antidot lattices as well as 
an indication of the range of validity of the less computationally expensive methods. 
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Figure 3. Band structures of three representative graphene antidot lattices. Full 
lines indicate results obtained by solving the Dirac equation (DE), while tight- binding 
results (TB) are shown with red dashed lines. Within these computational approaches 
we have exact particle-hole symmetry, and consequently only positive energies are 
shown. Note the different energy scale on the leftmost figure. 

Both the finite-element solutions of the Dirac equation (DE) and our tight-binding 
calculations (TB) were carried out on a standard PC, and a single band structure 
calculation could typically be performed in a few minutes for the relatively small-scale 
graphene antidot lattices considered in the following. The density functional theory 
calculations (DFT) were carried out on 8 AMD Opteron CPUs in parallel and typically 
lasted around 48 hours. Unlike the TB and the DFT methods, the computational time 
of our DE scheme does not increase with the size of the unit cell, determined by L, 
but only depends on the ratio R/L. For large unit cells, the DE scheme will therefore 
outperform both the TB and the DFT methods in terms of computational time. 

In Fig. [3] we show band structure results for three representative graphene antidot 
lattices using DE and TB. Both methods predict band gaps of a few hundred meVs 
for these relatively small dimensions of graphene antidot lattices. For low energies, DE 
predicts well the qualitative features of the bands obtained using TB, but the deviations 
become pronounced at higher energies. This is not surprising as the Dirac equation is 
only a valid description at low energies, where the band structure of bulk graphene is 
linear. Roughly, this means energies below 0.1/3 ~ 0.3 eV. Additionally, the increased 
kinetic energy due to the confinement of the particles renders the DE results less accurate 
for large antidot radii relative to the dimensions of the unit cell. This is apparent in 
the figure, where the bands at higher energies become increasingly inaccurate as the 
antidot radius is increased. However, even for the {L,R} = {10,6.4} structure, there 
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{12,3} 


{7,3} 


{10,6.4} 




eV 


A{i2,3} 


eV 


A{i2,3} 


eV 


A{12,3} 


DE 


0.54 (0.29) 


1 


1.27 (0.82) 


2.35 (2.83) 


1.53 (1.22) 


2.83 (4.21) 


TB 


0.23 


1 


0.74 


3.22 


1.01 


4.39 


DFT 


0.19 


1 


0.61 


3.21 


0.82 


4.32 



Table 1. Band gaps of three representative graphene antidot lattices. We show 
results obtained by solving the Dirac equation (DE), via tight-binding calculations 
(TB), and using density functional theory (DFT). Values in parentheses are obtained 
using DE and corrected for the low-radius behavior (see text). The band gaps are 
given in eV as well as in dimensionless values relative to the size of the band gap for 
the {L, R} = {12, 3} structure. 

is a qualitative agreement between the shapes of the bands at low energies found using 
the two methods. 

While the shapes of the bands at low energies are approximately the same for 
the DE and TB approaches, the sizes of the band gaps vary significantly. For the 
{L,R} = {12,3} structure, the band gap predicted by DE is more than twice as 
large as that obtained using TB. These differences may be traced back to two of the 
underlying assumptions of DE: linear bands of bulk graphene and absence of scattering 
between the two Dirac points. In order to illuminate the discrepancy we consider two 
limiting cases. We first consider the limit of large unit cells, i.e., large values of L. By 
investigating a large sample of different graphene antidot lattice using TB, some of us 
have demonstrated a simple scaling- law between the hole size and the band gap Eg, 
showing that Eg oc y/Ni^/Nceii for small values of the ratio R/L [22|. Here, A^hoie oc 
is the number of carbon atoms that have been removed from the intact unit cell in order 
to create the hole, and A^^ceii oc is the total number of carbon atoms in the intact unit 
cell (before the hole is made). We then find Eg oc {R/L)/L, showing that for a fixed 
value of the geometric ratio R/L, the band gap Eg falls off as 1/L with increasing L. 
For sufficiently large unit cells we thus expect the band gap to be well within the energy 
window for which the electronic bands of bulk graphene in fact are linear, and the band 
gaps obtained using DE should thus agree better with TB. The limit of small holes, i.e., 
R going to 0, is another important check point of our methods. In the DE approach, 
the boundary condition along the edge of a hole enforces a phase shift between the two 
spinor components, given by the angle (p indicated in Fig. [2b. With R going to 0, the 
phase shift must occur at a single point in space, resulting in a completely undetermined 
phase relationship at this point. Consequently, there is no adiabatic transition from a 
graphene antidot lattice to bulk graphene in the limit of vanishing hole sizes. Indeed, 
for small values of R, we find a non-vanishing band gap using DE, and extrapolating 
the results to i? = we find a band gap of the approximate size 1.02/3 /L. If we correct 
for this by simply subtracting this value from the band gaps calculated using DE, we 
find better agreement with the results obtained using TB, as shown in Table [H For 
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Figure 4. Band structures of three representative grapliene antidot lattices. Full lines 
indicate results obtained using density functional theory (DFT), while tight-binding 
results (TB) are shown with red dashed lines. Within the DFT scheme, particle-hole 
symmetry is not assumed, and we thus show results for energies both above and below 
the Fermi energy at zero. 



large values of L, the correction tends to zero, as we would expect. 

In Table [1] we also show results for the band gaps using density functional theory. 
The band gaps calculated using DFT are within 30% of the corresponding TB results, 
with DFT consistently reporting lower band gaps than TB. This follows the general 
tendency that energy gaps are underestimated in DFT [17] . For the three structures 
shown here, the band gaps increase with increasing relative hole size. This trend is 
captured well by all three methods. The band structures calculated with DFT are shown 
in Fig. m together with results obtained using TB for comparison. Generally, there is 
a reasonable qualitative agreement between the two methods in terms of the shapes of 
the bands, in particular, at energies close to the Fermi level. At larger energies, the 
qualitative features start to deviate significantly. Unlike the TB calculations, the DFT 
approach does not imply particle-hole symmetry, and the corresponding band structures 
are not symmetric around the Fermi level at zero. This difference is clearly seen in the 
figure. 

In contrast to the DE and TB calculations, our DFT scheme also includes the 
spin degree of freedom and is thereby able to predict the magnetic properties of the 
graphene antidot lattice. Within the TB description, graphene is considered a bipartite 
lattice structure with two sublattices, A and B, with non-zero hopping elements between 
different sublattice sites only. In that case, the total magnetic moment per unit cell M, 
can be determined from Lieb's theorem |48j, stating that M = Na—Nb with Na{b) being 
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Figure 5. Band structures of three representative graphene antidot lattices calculated 
with DFT. Full lines indicate results obtained using the DZP basis set, while dashed 
lines correspond to the smaller SZ basis set (see text). The real-space representations 
of the states corresponding to the points a and h are shown in Fig. [Sj 

the number of sites of sublattice A{B) in the unit cell. By inspection of the structures 
in Fig. [H we see that they have zero sublattice imbalance, i.e., Na = Nb, and we thus 
expect a zero total magnetic moment according to Lieb's theorem. Although our DFT 
calculations are not based on a description of graphene in terms of two sublattices with 
only nearest-neighbor coupling, we still find a zero total magnetic moment. Additionally, 
we find that no local magnetic moments are formed in any of the investigated structures. 
Lieb's theorem does not concern the formation of local magnetic moments, but the 
absence in the present cases can be understood from the circular shapes of the holes, 
which inhibit the formation of longer zig-zag shaped parts of the edge. This is similar 
to results obtained for graphene flakes, where relatively large zig-zag parts are needed 
for local magnetic moments to form [l9]. We thus conclude that the bands are all spin 
degenerate in the cases we have investigated. 

Within our DFT scheme, the computational time can be reduced by using a smaller 
basis set. We thus compare results obtained with the double-^ polarized (DZP) basis 
set involving 13 basis functions per carbon atom, used thus far, and results obtained 
using a single-^ (SZ) basis with only 4 basis functions per carbon atom. Results for 
the band structures obtained using the two different basis sets are shown in Fig. [5l 
The band structures obtained using the smaller SZ basis agree well with those obtained 
using the DZP basis, and the computational time is significantly reduced. An interesting 
difference between the band structures obtained using DFT compared to DE and TB, 
is the very low dispersion of the band roughly 0.5 eV below the Fermi level for the 
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Figure 6. Real-space representation of electronic states. The left panel corresponds 
to the point on the flat band in Fig. [5] indicated by a. For comparison, the right panel 
shows the state corresponding to the point b in Fig. \5\ We show the absolute square 
of the wavefunctions. 

{L,R} = {10,6.4} structure. The absolute square of the wavefunction for one of the 
spin degenerate states at the F-point, denoted by a in Fig. [5l is shown in Fig. [61 For 
comparison, we also show the state denoted by b in Fig. [H The state on the flat band 
is strongly localized on the zig-zag parts of the edge. The lower dispersion compared 
to TB is possibly due to the gradually increasing total electronic potential within DFT, 
when approaching the edge of a hole. The increased on-site energy of the edge atoms 
within DFT may thus cause stronger localization. 

5. Passivation 

Finally, we discuss the influence of edge passivation of the holes with hydrogen. In 
order to address this question we employ our DFT scheme. Details of the edges are not 
considered within our finite-element solutions of the Dirac equation (DE), and within 
a tight-binding description (TB), passivation is typically included simply as a shift of 
the hopping integral between carbon atoms along the edges due to the relaxed carbon- 
carbon bond length [50] . This correction leads to slightly increased energy gaps but has 
been ignored in the TB calculations in the present work. In contrast, DFT carefully 
treats the presence of hydrogen along the edge of a hole, and, importantly, the method 
includes the spin degrees of freedom, which turns out to be crucial in determining the 
influence of passivation on the electronic properties. We consider as an illustrative 
example the structure {L, R} = {4, 2} depicted in Fig. [3, shown with and without 
hydrogen passivation. We note that the hole geometry in this case is hexagonal, leading 
again to a vanishing magnetic moment without passivation. 

The resulting band structures, shown in Fig. [HI with and without passivation 
are very different. With full hydrogen passivation, several bands are spin degenerate. 
This degeneracy is lifted without passivation, and low dispersion bands stemming from 
dangling bonds are clearly present. The dispersions of these bands arise due to coupling 



Electronic properties of graphene antidot lattices 



15 



Figure 7. The unit cell of the {L, R} — {4, 2} structure. The structure is shown with 
(right) and without (left) complete hydrogen passivation of the carbon atoms along 
the edge of the hole. 



between neighboring edge atoms. Each danghng bond introduces a calculated spin of 
one Bohr magneton, l.OO/is, giving a total magnetization of 12.00/iB per cell, causing 
the lifting of the spin degeneracy. This magnetization involves only the sp^-orbitals 
and is strongly localized at the sites of the dangling bonds. We find that structural 
relaxation has no qualitative impact on the results in these two cases. 

Next, we investigate the effects of a single carbon vacancy at the edge. This 
introduces a sublattice imbalance of \Na — Nb\ = 1, resulting in an expected non-zero 
total magnetic moment according to Lieb's theorem. Elaborating on Lieb's work, Inui, 
Trugman and Abraham have shown that such a sublattice imbalance is accompanied 
by at least \Na — Nb\ midgap states with zero energy for a perfect bipartite lattice 
|51] . A recent discussion of similar statements can be found in Ref. [52]. In Fig. 
|9] we show the geometries of a single carbon vacancy at the edge, both with and 
without hydrogen passivation, as well as with and without having relaxed the geometries. 
The corresponding band structures are shown in Fig. [TOl Generally, we find two low 
dispersion bands close to the Fermi level. These are the midgap states with an induced 
spin splitting. The spin degeneracy is lifted for all bands due to the non-zero total 
magnetic moment. The main finding in the case without passivation of the atoms close 
to the vacancy are the two flat bands stemming from the dangling bonds, indicated 
with arrows in Fig. [TOl The dangling bonds are found to overlap in case of which it 
is energetically most favorable for the dangling bonds to have zero total spin as our 
calculations show. We find a magnetization of l.OOyUs per unit cell for both systems 
in the unrelaxed case. This magnetization is entirely due to the sublattice imbalance, 
and is, contrary to the case of dangling bonds, largely non-local, residing mainly on the 
TT-orbitals. 

Whereas relaxation has minor effects when passivation is included, the opposite 
is true for carbon vacancies without passivation. In that case, the two unpassivated 
carbon atoms at the edge approach each other, forming a pentagon as seen in Fig. 
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Figure 8. Band structure of the {L,R} — {4,2} graphene antidot lattice calculated 
with DFT. The left (right) panel shows the band structure without (with full) hydrogen 
passivation, corresponding to the unit cells in Fig. \7\ The systems are fully relaxed 
and the spin degree of freedom is included. Majority (minority) spin is shown with 
black (green). The Fermi level is at -E = 0. 



[9li. A similar phenomenon has been observed theoretically for single carbon vacancies 
in bulk graphene, where the spin of such vacancies can usually be understood as an 
unsaturated dangling bond on the neighboring carbon atom, not forming the pentagon 
[53| [Ml [55] . This results in a calculated magnetic moment of around l/i^. In our case, 
however, there are only two neighboring atoms. In fact, in both cases the magnetic 
moment is better understood using Lieb's theorem, as discussed by Palacios, Fernandez- 
Rossier, and Brey, in the case of carbon vacancies in bulk graphene [56]. In the pentagon 
geometry, two sites belonging to the same sublattice bond stronger to each other, which 
is reflected in the smaller bond length of 1.67 A compared to 2.46 A in the case without 
relaxation. Consequently, the dangling bonds are then saturated and the corresponding 
flat bands are not present. Additionally, the bipartite lattice symmetry is broken, 
causing a reduction in the magnetic moment from l.OO/is to 0.50fiB- Consequently, 
the spin splitting of the bands is reduced. The midgap states are still observed, but in 
this case with more dispersion. The features of the bipartite lattice are thus maintained 
in a moderated version, when pentagons are formed due to a carbon vacancy along the 
edge of the hole. We stress that while the magnetization arising from Lieb's theorem 
is predictable, the magnetization due to dangling bonds is merely a result of energy 
optimization and therefore harder to predict. 
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Figure 9. Single carbon vacancy at tiie edge of the liole in the {L, i?} = {4, 2} 
structure. In the left (right) panels the structures have not been (have been) relaxed. 
In the upper (lower) panels, the carbon atoms next to the vacancy have (have not) 
been passivated with hydrogen. The calculated magnetic moment is indicated in each 
panel. 

6. Conclusions 

We have carried out a numerical study of the band structures of graphene antidot 
lattices, using three different computational approaches of varying levels of complexity 
and accuracy. Finite-element solutions of the Dirac equation (DE) provide a simple and 
fast scheme, capturing essential qualitative features of the band structures and band 
gaps. For more reliable predictions of the band structures, we employed a 7r-orbital 
tight-binding scheme (TB) as well as computationally heavy density functional theory 
calculations (DFT). The three methods all predict an opening of a band gap on the order 
of a few hundred meVs for the nano-scale structured graphene antidot lattices studied 
in this work. Qualitative similarities were found for the band structures calculated with 
the three different methods. Finally, we discussed the effects of hydrogen passivation 
along the edges of the holes. Passivation was found to have a significant influence on the 
band structures, and the presence of carbon vacancies along the hole edges were shown 
to induce midgap bands. 
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Figure 10. Band structures of the {L, R} — {4, 2} graphene antidot lattice with 
a single carbon vacancy in the unit cell. Dashed lines indicate band structures for 
the unrelaxed geometry shown in Fig. [9l panels (a) and (c), while full lines are the 
corresponding results for the relaxed structures, panels (b) and (d). The unfilled 
bands of the dangling bonds are indicated in the left panel by horizontal arrows. 
The corresponding filled bands at lower energies are not shown in the plot. Majority 
(minority) spin is shown with black (green). The Fermi level is at E ~ 0. 
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